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Passive  localization  with  propagation 
models 

Edmund  J.  Sullivan 


Abstract:  Passive  ranging  using  propagation  models  can  be  effected  in 
two  different  ways:  forward  modeling  and  direct  inversion.  The  forward¬ 
modeling  method  performs  an  exhaustive  recomputation  of  the  field  over  a 
set  of  assumed  source  locations  and  seeks  the  best  match  to  the  measured 
data.  The  direct  inversion  method  inverts  a  set  of  linear  equations  consti¬ 
tuted  by  the  normal  mode  model.  An  overview  of  these  two  approaches  is 
given  along  with  a  summary  of  studies  performed  to  date.  This  is  followed 
by  the  presentation  of  a  maximum-likelihood  algorithm  for  the  direct  in¬ 
version  method.  A  discussion  of  the  merits  and  weaknesses  of  each  method 
is  given.  Also,  an  appendix  which  discusses  some  general  aspects  of  the 
inverse  problem  is  included. 


Keywords:  bearings-only  o  deep  water  o  environmental  acoustics  o 
multipath  ranging  o  passive  localization  o  propagation  models  o  sonar 
o  target  tracking  o  towed  array  o  wavefront  curvature 
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1.  Introduction 


One  means  of  obtaining  improved  processing  performance  for  sonar  systems  is  to 
include  environmental  information  in  the  processing  scheme  [1-3].  In  spite  of  the 
potential  gains  that  this  would  offer  however,  it  has  only  recently  begun  to  be 
exploited.  There  are  good  reasons  for  this.  One  is  that  in  most  cases  of  inter¬ 
est  sufficiently  accurate  environmental  information  has  simply  not  been  available. 
A  second  reason  is  that  the  processing  load  demanded  by  such  a  scheme  can  be 
prohibitive.  One  area  where  these  obstacles  are  slowly  being  overcome  is  that  of 
passive  localization. 

The  standard  localization  techniques,  which  usually  concentrate  on  range  and  bear¬ 
ing,  either  do  not  use  environmental  information  or  include  it  in  a  minimal  sense. 
These  standard  techniques  fall  into  three  general  categories:  bearings-only,  wave- 
front  curvature  and  multipath. 

The  bearings-only  approach  consists  of  making  at  least  two  independent  bearing 
estimates  with  a  towed  array.  The  range  is  then  computed  from  these  two  estimates 
[4].  More  general  target  tracking  schemes  based  on  towed  array  bearing  estimates 
have  been  developed  which  provide  more  complete  target  motion  analyses.  These 
techniques,  however,  can  be  very  costly  in  terms  of  processing  time.  The  wavefront 
curvature  technique  provides  an  estimate  of  range  as  the  radius  of  curvature  of  the 
wavefront  under  the  assumption  that  the  wave  is  cylindrical  [5].  This  can  be  done 
either  with  a  towed  array  or  at  least  three  independent  hydrophones  in  a  horizontal 
line.  The  quality  of  the  estimate  depends  strongly  on  the  size  of  the  horizontal 
acoustic  aperture  at  right  angles  to  the  direction  of  propagation.  Both  the  bearings- 
only  and  the  wavefront  curvature  techiques  are  severely  limited  by  the  horizontal 
acoustic  aperture  available  whereas  the  third  technique,  that  of  multipath  ranging, 
depends  on  vertical  aperture  or  horizontal  aperture  along  the  bearing  line,  since 
it  relies  on  the  ability  to  resolve  multipath  components  of  the  fields  [6].  It  should 
be  pointed  out  that  this  technique  does  depend  on  environmental  information.  In 
particular,  one  must  know  the  water  depth  and  the  sound  velocity  profile  (SVP). 
Thus,  this  technique  uses  a  deep-water  propagation  model  to  provide  an  estimate 
of  range.  Like  the  bearings-only  technique,  it  is  based  on  a  plane-wave  assumption. 
The  wavefront  curvature  technique  uses  the  non-planar  character  of  the  wave  under 
the  assumption  that  it  is  a  circular  wave  but  uses  no  environmental  information. 

The  techniques  that  we  wish  to  discuss  in  this  report  constitute  the  next  evolution¬ 
ary  step  in  that  they  include  the  environmental  information  in  a  more  complete 
and  more  consistent  manner.  They  also  differ  from  the  above  techniques  in  that 
the  studies  carried  out  up  to  the  present  are  based  on  shallow- water  models. 
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2.  Theory 


2.1.  BACKGROUND 

If  one  has  a  model  that  describes  the  acoustic  propagation  in  a  given  channel,  then 
given  sufficient  information  about  an  acoustic  source  in  this  channel  and  enough 
sufficiently  accurate  environmental  information  to  provide  the  proper  parameters 
for  the  model,  the  acoustic  field  can  be  computed.  The  accuracy  of  the  field  as 
computed  from  the  model  depends  on  two  things:  the  goodness  of  the  environmental 
parameters  and  the  fidelity  of  the  model.  This  basically  describes  what  is  known  as 
the  ‘forward’  problem.  That  is,  given  sufficient  information  about  the  source  and 
the  medium,  compute  the  resulting  field.  The  localization  techniques  constitute  an 
‘inverse’  problem:  Given  sufficient  information  concerning  the  acoustic  field  and  the 
medium,  compute  the  position  of  the  source.  This  we  refer  to  as  an  inverse  problem 
of  the  first  kind.  The  inverse  problem  of  the  second  kind  obtains  when  the  field  and 
the  source  are  known  and  the  properties  of  the  medium  are  to  be  computed.  This 
would  be  of  interest  in  the  case  of  active  sonar,  where  one  would  be  concerned  with 
determining  certain  characteristics  of  the  channel.  Here,  we  are  concerned  with  the 
inverse  problem  of  the  first  kind.  The  inverse  problem  entails  certain  mathematical 
difficulties  that  do  not  appear  in  the  forward  problem.  First,  the  solution  may  not 
be  stable  to  small  changes  in  either  the  data,  the  model  parameters  or  the  noise. 
Also,  the  solution  may  not  be  unique.  If  either  or  both  of  these  situations  occur, 
the  problem  is  said  to  be  ill-posed  [7].  (A  more  complete  discussion  of  the  inverse 
problem  can  be  found  in  the  appendix.)  The  importance  of  this  to  the  passive 
localization  problem  is  that  for  each  algorithm  and  model  under  consideration,  a 
careful  sensitivity  study  should  be  carried  out  in  order  to  determine  the  effects  of 
measurement  and  parameter  errors  on  the  localization  estimate. 

The  present  efforts  in  model-based  passive  localization  are  focussed  on  two  ap¬ 
proaches:  matched-field  processing  and  direct  inversion.  These  will  be  discussed  in 
turn,  beginning  with  matched-field  processing. 


2.2.  MATCHED  FIELD  PROCESSING 


Matched-field  processing  is  best  described  as  the  solution  to  an  inverse  problem 
by  forward  modeling.  Suppose  the  range  and  depth  of  a  point  acoustic  source  are 
to  be  estimated  from  the  data  received  on  a  vertical  array  of  hydrophones.  The 
field  as  predicted  by  the  model  for  a  source  location  at  some  arbitrary  point  on 
a  specified  range-depth  grid  is  computed.  This  is  then  compared  to  the  measured 
field  in  some  manner.  For  example,  the  estimator  could  be  the  correlation  between 
the  two  fields.  This  estimator  is  then  computed  for  all  range-depth  combinations 
and  plotted  on  a  range  depth  map.  The  location  of  the  best  estimate  on  this  map 
then  constitutes  the  estimate  of  the  true  range  and  depth  of  the  source. 
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An  example  of  such  an  estimator  is  Bucker’s  detection  factor  [8].  This  can  be 
written  as 

D({a})  =  Pl({a})RPM(  {a}),  (1) 

where  R  is  the  covariance  matrix  of  the  measured  field  and  Pm  ({a})  is  the  field 
as  predicted  by  the  model  for  source  coordinates  {a}.  For  a  vertical  array  of  N 
phones,  Pm({ol})  is  an  N- dimensional  complex  vector  and  R  is  N  x  N.  (It  should 
be  noted  here  that  in  Bucker’s  original  definition  the  diagonal  of  R  is  removed.) 
Equation  (1)  can  be  thought  of  as  the  power  output  of  a  beamformer  ‘steering 
vector’  PM({a})- 

The  detection  factor  is  maximum  when  the  modeled  field  equals  the  measured  field. 
This  can  be  seen  as  follows.  Assuming  that  the  estimate  of  the  covariance  matrix 
is  obtained  by  a  time  average,  we  have 

D  =  PlRPM  =  (2) 

where  we  have  temporarily  dropped  the  explicit  dependence  on  {a}  for  clarity. 
Since  Pm  is  deterministic  (for  the  case  of  a  deterministic  model),  Eq.  (2)  can  be 
written  as 

D  =  (PmPP'Pm)  =  ((Pl,P)(Pl,PV).  (3) 

By  the  Schwartz  inequality,  the  inner  product  Pj^P  is  maximum  when  Pm  —  P. 

More  recent  studies  have  considered  other  forms  for  the  detection  factor  or  ‘ambi¬ 
guity  function’  as  it  is  called  by  many  authors.  These  efforts  are  based  on  replacing 
the  ‘beamformer’,  Pm ,  in  Eq.  (1)  with  a  high-resolution  beamformer.  In  particular, 
attention  seems  to  be  focussing  on  the  so-called  maximum-likelihood  beamformer 
of  Capon  [9].  This  is  a  beamformer  whose  weights  are  such  that  the  spatial  spec¬ 
trum  is  a  maximum- likelihood  estimate  of  the  power  being  received  as  a  function  of 
wavenumber  (angle)  if  the  noise-only  covariance  matrix  is  known  [10].  It  is  equiva¬ 
lent  to  the  beamformer  that  obtains  if  the  power  is  minimized  everywhere  except  in 
the  look  direction,  which  is  constrained  to  be  a  constant.  The  resulting  beamformer 
weights  are  given  by 


R-1Pm({*}) 

Pl({a})R-ipM({a})' 


Substituting  Eq.  (4)  into  Eq.  (1)  yields  for  the  ambiguity  function 


(4) 


1 


(5) 


Several  investigators  have  used  this  form  to  estimate  range  and  depth  [11-14].  The 
results  indicate  that  the  lower  sidelobes  produced  by  this  technique  tend  to  reduce 
the  aliasing  problem  at  the  cost  of  some  detection  performance. 
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There  tends  to  be  some  confusion  in  the  hterature  as  to  just  what  Eq.  (5)  con¬ 
stitutes.  That  is,  since  it  is  based  on  a  maximum- likelihood  beamformer,  does  it 
constitute  a  maximum- likelihood  estimate  of  the  source  position?  What  Eq.  (5) 
does  constitute  is  a  maximum- likelihood  estimate  of  the  beamformer  output  power 
for  any  parameter  set  {a}.  The  matched-field  technique  then,  is  a  graphic  or  numer¬ 
ical  technique  which  estimates  {a}  by  finding  the  maximum  of  this  power  estimate. 
It  is  by  no  means  clear,  however,  that  this  criterion  delivers  a  maximum  likelihood 
estimate  of  source  position. 

2.3.  DIRECT  INVERSION  TECHNIQUE 


The  direct  inversion  technique  is  based  on  the  fact  that  the  normal  mode  model  of 
propagation  permits  a  set  of  linear  equations  to  be  written  that  can  be  inverted. 
The  source  position  coordinates  are  then  extracted,  usually  as  a  second  step.  The 
technique  is  also  sometimes  referred  to  as  a  modal  filter. 


The  normal  mode  model  is  expressed  as  follows: 

N 


p(v  y\  __  ,  .2  2  1  <ftn(z)<ftn(^0 )  ~anr  +  tfcnr-i7r /4 

r ,  z)  —  w  P  0  rz 7,  n —  e  ’ 

Vo7rr  VK 


(6) 


where  P(r,  z)  is  the  acoustic  pressure  at  depth  2  and  range  r  from  a  point  source  at 
depth  zo,  u>  is  the  circular  frequency,  po  is  the  water  density,  an  is  the  loss  factor  for 
mode  n,  and  kn  is  the  horizontal  wavenumber  which  satisfies  u >2/c2  =  fc2 +72.  Here, 
c  is  the  speed  of  sound  and  can  depend  on  depth.  7n  is  the  vertical  wavenumber. 
Equation  (6)  is  the  solution  by  separation  of  variables  to  the  Hemholtz  form  of  the 
wave  equation.  Using  a  far-field  assumption  in  shallow  water,  the  field  is  considered 
to  be  a  cylindrical  wave.  Thus  the  total  field  is  considered  to  be  a  function  of  r  and 
2  only.  Upon  separation,  the  equation  in  z  becomes  an  eigenvalue  equation  given 

by 


d2<f>n 

dz 2 


(7) 


whose  eigensolutions  are  the  so-called  modal  functions.  Upon  solving  the  equa¬ 
tion  in  r,  which  yields  a  Hankel  function  as  the  range-dependent  part  of  the  total 
solution,  and  taking  the  far-field  approximation  to  the  Hankel  function,  Eq.  (6) 
obtains. 


A  more  comprehensive  discussion  of  the  normal-mode  propagation  model  can  be 
found  in  [15]. 


For  our  purposes,  knowledge  of  the  absolute  value  of  the  field  is  not  necessary,  thus 
we  normalize  Eq.  (6),  yielding 


P(r,z) 


N 


E 


<Pn(z)<t>n{zo)  -anr  +  iknr 

VK 


(8) 
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Suppose  now  that  we  wish  to  determine  the  position  coordinates  of  the  point  acous¬ 
tic  source  by  means  of  a  vertical  array  of  hydrophones.  Defining 


Eq.  (8)  becomes 


m  ,n  —  4*n{Zm)  /  \/^nj 

(9a) 

An  -  <^n(Zo)e-^r, 

(9b) 

v  —  4  J^rtT 

An  —  > 

(9c) 

Pm  —  Af  m ,  n  X  n  ♦ 

(10) 

Inversion  of  Eq.  (10)  yields  the  N  functions  \n  which  in  principle  can  be  solved 
for  the  range  and  depth,  since  Mm>n  and  kn  (and  also  an)  can  be  determined  from 
knowledge  of  the  depth  of  the  water,  the  bottom  conditions  and  the  sound  velocity 
profile;  that  is,  without  knowledge  of  the  source  coordinates.  The  problem  can  be 
generalized  somewhat  by  the  addition  of  noise.  Defining  a  noise  vector  £n,  we  have 


P m  —  m ,  n  X  n  £  m  >  1  ^  ^  ^  ^  i 

1  <m<M,  (11) 


where  M  is  the  number  of  hydrophones  and  N  is  the  number  of  modes  supported 
by  the  acoustic  channel. 

Solving  Eq.  (11)  by  the  method  of  maximum- likelihood,  the  maximum- likelihood 
estimation  of  Xn  is,  for  the  case  of  M  >  N ,  given  by  [16] 

x  =  {M' R~l  M)~l  M' R~l  P,  (12) 


where  we  have  dropped  the  matrix  and  vector  indices  for  convenience.  Equa¬ 
tion  (12)  is  valid  for  gaussian  noise  and  is  sometimes  referred  to  as  the  generalized 
least-squares  estimator,  R  is  the  noise  covariance  matrix  and  is  given  by 

R=  E{ee ♦}.  (13) 


Hinich  used  this  approach  to  develop  a  technique  for  the  estimation  of  the  depth 
of  a  point  source  with  a  vertical  array  [17].  This  is  probably  the  first  time  that  a 
sophisticated  propagation  model  was  used  for  source  localization. 

More  recently,  Shang,  Clay  and  Wang  developed  a  technique  for  range  estimation 
based  on  Eq.  (10)  [18].  Their  work  is  based  on  the  observation  that,  since  the  modal 
matrix  Mm>n  is  real,  it  is  only  necessary  to  know  the  phase  of  the  modal  functions 
Xn  to  determine  range.  Thus,  Eq.  (10)  is  solved  directly  for  the  functions  \n.  From 
Eq.  (9c),  one  then  has 

arg(Xn)  -  arg(Xm)  =  $n,m,  (14) 
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where  the  modal  phase  difference  $U|m  is  given  by 

^n,m  —  ( kn  ^m)**  Lnrn‘2'K.  (15) 

The  term  LntTn arises  since  $n>m  only  contains  information  concerning  the  prin¬ 
cipal  values  of  the  phases. 

A  third  mode  is  then  introduced  permitting  the  calculation  of  a  second  modal  phase 
difference,  say  where 


^n,p  —  i^n  Lrxp<l'K. 


(16) 


The  range  is  then  eliminated  between  Eqs.  (15)  and  (16)  resulting  in 


^n,m  ,  j  (  kn  ~  km  \ 

27T  +  n>m  "Un-tp/ 


(17) 


A  search  is  then  made  to  find  all  integers  £n>m  and  LntP  that  satisfy  the  condition 
that  the  integer  parts  and  the  fractional  parts  of  Eq.  (17)  are  simultaneously  equal. 
All  of  the  solutions  of  Eq.  (17)  are  not  valid  since  the  sign  of  <J>nm  is  not  known. 
Shang  et  al.  discuss  a  method  to  deal  with  this  problem. 


A  natural  generalization  of  the  work  of  Shang  et  al.  is  effected  as  follows.  Starting 
from  Eq.  (12),  a  solution  can  be  formally  written  as 


x  =  VP0  +  Ve  =  VP,  (18) 

where  Pq  is  the  noise-free  signal  and 

V  =  {M*R-lM)~lM*R-1.  (19) 

Equation  (19)  is  taken  from  Eq.  (12).  That  is,  Eq.  (18)  is  the  generalized  least- 
squares  solution  to  the  problem.  A  ‘solution  matrix’  is  then  defined  by 


s  =  (xx1).  (20) 

where  the  brackets  indicate  the  time  average.  For  the  case  of  noise  uncorrelated 
with  the  signal,  substitution  of  Eq.  (18)  into  Eq.  (20)  yields 

5  =  V(P0P£)V'  +  VRV\  (21) 

where  we  have  used  Eq.  (13). 

This  formulation  has  several  advantages  associated  with  it.  First,  it  allows  the 
data  to  be  introduced  in  terms  of  the  covariance  matrix,  which  eliminates  the  time 
dependence  from  the  problem.  Second,  it  permits  numerical  studies  of  various  types 
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of  noise  to  be  directly  carried  out  in  terms  of  the  noise  covariance  matrix.  Third, 
the  arguments  of  the  elements  of  the  solution  matrix  directly  provide  the  phase 
differences  necessary  for  the  solution. 

Numerical  studies  of  the  algorithm  of  Shang  et  al.  (Eq.  (IT))  indicate  that  it  is  quite 
sensitive  to  errors  in  the  positions  of  the  hydrophones.  A  more  stable  solution  can 
be  found  as  follows.  From  the  O'th  element  of  the  solution  matrix  S,  one  can  write 

arg(Sfj)  =  (ki  —  kj)r  -  Nir.  (22) 

The  term  N tv  arises  from  the  fact  that  arg(5;_,)  is  only  the  principal  value  of  the 
phase,  and  the  true  sign  of  Xn  is  not  known. 

Solving  for  r,  the  result  is 


r 


ij 


arg  {Sjj) 

{ki  -  kj) 


+  N 


{ki  —  kj) 


With  the  definitions 


o  arg  {Sjj) 

ij  (ki-kjY 


A 


( ki-kj ) 


Eq.  (23)  becomes 

rH  =  r°ij  +  NArtj. 


(23) 


(24) 

(25) 


By  introducing  a  third  mode,  two  more  solutions  for  the  range  can  be  written,  say 
rih  and  rjfc.  If  now  all  the  multiple  solutions  for  each  of  the  three  range  solutions 
are  computed  out  to  some  maximum  range,  these  solutions  can  be  counted  in  range 
bins.  At  the  bin  containing  the  correct  range  one  would  expect  three  solutions. 
This  can  be  easily  generalized  to  as  many  modes  as  described.  For  example  using 
five  modes,  10  pairs  can  be  constructed  so  that,  at  the  range  bin  containing  the 
correct  range  10  solutions  should  appear.  A  histogram  can  then  be  constructed 
indicating  the  number  of  solutions  in  each  range  bin.  The  range  corresponding 
to  the  bin  containing  the  largest  number  of  solutions  is  the  estimated  solution. 
Numerical  studies  have  indicated  that  this  algorithm  is  reasonably  stable  to  errors 
in  the  hydrophone  positions. 

Once  the  range  is  determined,  the  depth  can  be  found  by  using  Eqs.  (9b),  (9c) 
and  (12).  The  procedure  is  as  follows.  Equation  (9b)  is  used  to  eliminate  An  from 
Eq.  (9c).  The  result  is 

x  =  Mzo)e{kn+ian)f  ■  (26) 

Here,  \n  is  the  estimate  of  \n  taken  from  Eq.  (12).  The  complex  wavenumbers 
kn  +  ian  are  known  and  f  is  the  range  estimate.  Thus,  the  only  unknown  in 
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Eq.  (26)  is  z0 ,  the  source  depth.  It  can  be  seen  that  the  solution  of  Eq.  (26) 
requires  samples  of  the  modal  functions  for  all  values  of  the  desired  depth  mesh. 
However,  the  evaluation  of  the  matrix  M  (see  Eq.  (9a))  already  provides  a  set  of 
samples  at  the  hydrophone  depths.  Moreover,  the  SNAP  program  [19],  as  a  matter 
of  course,  computes  the  modal  functions  for  up  to  501  depths  uniformly  distributed 
over  the  water  column. 

A  solution  to  Eq.  (26)  by  the  method  of  least  squares  takes  the  following  form: 

N 

-  <M2m)]2  =  Jm-  (27) 

n=l 

The  quantity  Jm  can  be  plotted  as  a  function  of  m  with  its  minimum  providing 
the  estimate  of  source  depth.  The  vector  <£n  is  the  estimate  of  ^n(z0)  taken  from 
Eq.  (26)  and  </>n(zm)  is  the  matrix  of  computed  values  of  the  modal  functions. 
N  is  the  number  of  modes.  In  practice,  the  quantities  </>n  and  <?!>n(zm)  must  be 
normalized,  since  the  absolute  value  of  \  is  unknown. 


3.  Experimental  work 


There  have  been  several  successful  full-scale  localization  experiments  to  date,  all 
using  the  matched-field  technique.  As  for  the  direct  inversion  method,  there  has 
been  one  successful  scale-model  experiment. 

Fizell  and  Wales  [12]  report  a  successful  localization  using  a  980  m  array  in  the 
upper  part  of  the  water  column  in  water  of  depth  somewhat  greater  than  2800  m. 
This  experiment  took  place  in  an  arctic  environment,  where  both  the  array  and 
the  source  were  suspended  from  the  ice  cover.  The  20  Hz  source  was  correctly 
located  at  a  range  of  270  km.  The  signal  to  noise  ratio  was  quoted  as  ‘high’.  The 
propagation  model  used  was  a  fast  field  program. 

Feuillade  and  Kinney  [13]  report  on  an  experiment  performed  off  Panama  City, 
Florida,  involving  a  450  Hz  source  at  a  range  of  2.2  km  and  a  depth  of  31  m  in 
31.5  m  of  water.  The  array  was  approximately  26  m  in  length.  The  model  used 
was  SNAP  [19].  The  results  indicated  a  biased  range  estimate  (2.42  km)  and  some 
ambiguity  due  to  multiple  peaks  of  the  ambiguity  surface. 

Bucker  [20]  reports  a  successful  localization  experiment  using  a  3  kHz  source  in 
water  with  a  depth  of  66  ft.  He  reports  good  localization  out  to  ranges  of  at  least 
four  times  the  water  depth. 
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Finally,  Shang  et  al.  [18]  report  on  a  scale-model  experiment  with  a  3500  Hz  source 
at  a  range  of  1.8  m  in  air.  The  acoustic  channel  was  a  waveguide  formed  by  two 
parallel  plaster  plates  with  a  separation  of  0.123  m.  The  experimentally  determined 
range,  using  their  direct  inversion  approach,  was  1.798  in.  Three  modes  were  used 
in  the  calculation.  The  vertical  array  was  actually  a  single  receiver  which  was  used 
as  a  synthetic  array. 


4.  D  iscussion 

As  pointed  out  in  the  introduction,  from  a  logical  point  of  view,  both  the  matched- 
field  processing  method  and  the  direct  inversion  method  are  considered  inverse 
problems.  However,  in  practice  only  the  direct  inversion  method  is  considered  as 
an  inverse  problem.  The  reason  for  this  is  that  much  of  the  literature  concerning 
inverse  problems  deals  with  the  mathematical  issues  that  are  unique  to  the  method 
of  direct  inversion.  These  issues  are,  to  a  large  degree,  concerned  with  the  so-called 
‘ill-posed’  problem.  A  more  complete  discussion  of  this  is  given  in  the  appendix.  It 
is  of  concern  to  us  here  since  it  manifests  itself  in  the  method  of  direct  inversion. 
From  Eq.  (12),  we  see  that  the  solution  requires  the  inversion  of  a  matrix.  Even  in 
the  case  of  white  noise  where  the  covariance  matrix  drops  out  we  must  invert  the 
matrix  M*  Af .  When  the  data  do  not  provide  a  complete  enough  sample,  this  matrix 
becomes  singular.  In  particular,  when  the  array  is  vertical  studies  have  shown  that, 
at  least  at  low  frequencies,  nearly  the  complete  water  column  must  be  sampled  to 
stabilize  the  problem  [21,  22].  In  the  matched-field  technique  this  does  not  occur.  A 
more  ‘graceful  degradation’  of  the  precision  of  the  estimate  occurs  instead.  Another 
advantage  of  the  matched-field  method  is  that  one  is  not  limited  to  the  normal¬ 
mode  model.  In  principle  any  model  of  any  degree  of  sophistication  can  be  used.  In 
particular,  Baggeroer  et  al.  [11]  use  a  full-field  model  which  includes  the  continuous 
part  of  the  spectrum  as  well  as  the  modal  solution  [23],  The  direct  inversion  method, 
on  the  other  hand,  although  limited  to  the  range  dependent  case,  allows  a  direct 
solution  to  be  computed  without  an  exhaustive  recomputation  of  the  field.  In 
particular,  when  only  the  range  is  desired,  the  computation  consists  of  a  single 
matrix  inversion  (Eq.  (12))  followed  by  the  evaluation  of  Eq.  (23)  for  all  relevant 
values  of  N .  The  determination  of  depth  however  requires  a  search  over  source 
depth.  Thus  a  two-dimensional  search  is  replaced  by  two  one- dimensional  searches. 
Another  possible  advantage  of  the  direct  inversion  method  is  its  performance  against 
noise.  The  matched-field  techniques  appear  to  require  S/N  values  as  high  as  10  dB 
where  a  study  of  the  S/N  performance  of  the  direct  inversion  method  produced 
solutions  with  S/N  ratios  as  low  as  —13  dB  [22]. 

As  pointed  out  in  Subsect.  2.3,  the  solution  of  Eq.  (12)  provides  a  maximum- 
likelihood  solution  to  the  problem  for  the  case  of  gaussian  noise.  This  can  be  seen 
as  follows.  Equation  (12)  constitutes  a  linear  problem  in  the  vector  where  x 
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is  given  by  Xn  =  Anelknr  (Eq.  (9c)).  Thus,  solution  of  Eq.  (12)  provides  a  direct 
maximum-likelihood  estimation  of  x,  which  we  designate  as  x.  Maximum  likelihood 
estimators  possess  the  property  of  invariance  under  parameter  transformation  [24]. 
This  means  that  if  a  is  the  parameter  of  interest  and  we  have  the  maximum- 
likelihood  estimate  of  x,  designated  as  x,  then 


X  =  X(«)-  (28) 

This  can  be  seen  from  the  following  argument.  Given  the  appropriate  likelihood 
function  L(a)  (or  its  logarithm)  the  estimate  of  a  is  the  solution  of 


dL{a) 

da 


=  0. 


(29) 


But  this  can  be  written  as 


dL_  _  frLdZ 
da  ~  dZ  da' 


(30) 


where  Z(a)  is  a  parameter  that  permits  a  more  convenient  solution.  As  can  be  seen 
by  inspection  of  Eq.  (30),  that  value  of  a  that  satisfies 


(31) 


automatically  satisfies  Eq.  (30).  Thus, 


Z  =  Z(a ),  (32) 

which  is  equivalent  to  Eq.  (28).  Of  course,  the  question  of  whether  one  has  a  global 
maximum  or  not  still  remains. 
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Appendix  A 
Inverse  problems 

Generally,  acoustic  signal  processing  is  concerned  with  both  detection  and  estima¬ 
tion.  Detection  is  the  determination  of  the  existence  or  non-existence  of  a  signal 
at  the  receiver  whereas  estimation  is  the  determination  of  certain  parameters  or 
descriptors  of  the  signal,  the  medium  or  its  contents.  Certain  estimation  problems 
fall  into  a  class  of  problems  commonly  referred  to  as  'inverse  problems’.  This  is  not 
a  well-defined  concept  and  is  best  defined  in  terms  of  its  related  ‘forward  problem’. 
The  problem  under  discussion  here  is  the  passive  ranging  problem.  The  forward 
problem  then,  would  consist  of  computing  the  acoustic  field  radiated  from  a  source 
of  known  location  in  a  known  acoustic  channel.  The  concomitant  inverse  problem 
in  one  form  is  concerned  with  the  determination  of  the  location  of  the  source,  based 
on  the  measurements  of  the  acoustic  field  and  the  known  properties  of  the  acoustic 
channel.  The  inverse  problem  will,  almost  always,  entail  much  more  computational 
labour  than  the  forward  problem.  Most  linear  inverse  problems  can  be  written  in 
the  form  of  a  Fredholm  equation  of  the  first  kind: 

f(X)  =  J  KiX^'MX^dX',  (Al) 

or  in  its  discrete  form 

fi  =  Kij9j.  (A2) 

If  Eq,  (Al)  represented  an  acoustic  scattering  problem,  then  K(X,  X1)  would  rep¬ 
resent  the  acoustic  energy  scattered  from  a  differential  volume  dX1  or  ‘source’  at 
the  point  denoted  by  the  vector  coordinate  X1  and  g(X)  would  represent  the  dis¬ 
tribution  of  sources  and  therefore  the  distribution  of  scattering  medium.  The  scat¬ 
tered  field  would  be  given  by  the  function  f(X ).  The  forward  problem,  then, 
would  be  concerned  with  determining  f(X)  given  some  information  about  g(X) 
and  K(X1XI).  An  example  of  an  inverse  acoustic  scattering  problem  would  be 
the  determination  of  <7(A)  given  some  information  about  K{X)  X1)  and  measured 
values  of  f(X).  Sometimes  K(X)XI)  is  called  Green’s  function  and  here  would 
describe  the  acoustic  field  scattered  from  a  point  ‘source’. 

In  the  case  of  the  passive  localization  problems,  Eq;  (21)  reduces  to  the  form  of 
Eq.  (10)  where  one  makes  the  identification 

fi  =  Pu  (A3) 

Mij  =  Kij,  (A4) 

9i  =  Xi •  (A5) 

Because  of  its  generality,  the  inverse  problem  does  not  lend  itself  to  a  formal  struc¬ 
ture  in  the  same  sense  as  do  detection  and  estimation.  It  is  probably  more  realistic 
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to  consider  it  as  a  collection  of  techniques.  There  are,  however,  some  insights  that 
can  be  gleaned  from  examination  of  Eq.  (Al).  First,  it  is  clear  that  it  is  generally 
an  integral  equation  problem.  Except  for  rare  special  cases,  solution  in  closed  form 
will  not  be  possible,  so  that  a  numerical  solution  will  be  called  for.  If  the  data  are 
in  discrete  form,  it  becomes  a  problem  in  linear  equations.  Also,  there  can  arise 
severe  mathematical  pathologies  in  the  solution  of  such  problems.  One  of  these 
is  that  the  problem  can  be  what  mathematicians  refer  to  as  ill-posed.  This  is  a 
problem  that  can  be  due  to  noisy  data,  measurement  errors,  not  enough  data,  or  a 
combination  of  these.  A  problem  is  well  posed  when  a  unique  solution  exists  that 
is  stable  to  changes  in  the  data.  The  problem  of  nonuniqueness  can,  in  many  cases, 
be  dealt  with  by  introducing  a  priori  information. 

The  stability  problem  arises  when,  given  the  existence  of  a  solution,  the  solution  is 
extremely  sensitive  to  small  perturbations  in  the  data.  As  an  example,  consider  a 
measurement  of  f(X)  denoted  by  f(X i).  Let  this  measurement  incur  an  error  A / 
and  let  /0  be  the  true  value.  Then, 


Subtracting  Eq.  (A7)  from  Eq.  (A6)  yields 


(AS) 


It  can  be  seen  from  Eq.  (A8)  that  A /  can  be  thought  of  as  a  weighted  average  of 
A g  where  K(X\,  Ar)  is  the  weight.  Thus,  one  is  free  to  select  a  function  A g  whose 
weighted  average  is  as  close  to  zero  as  desired,  but  can  produce  large  errors  on 
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